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ABSTRACT 


A sea, swell and surf program is improved, tested and evaluated on a 
micro-computer (HP-9845B). Sea swell is ocr by a two dimensional 
spectral model. The energy balance equation is tested for different cases of 
wind velocities and water depths. Satisfactory agreement is observed 
between the offshore model and expected wave heights for a 15 knot wind, 
but the model overbuilds wave energy for a 30 knot’ wind. Wave 
transformation is described by a one dimensional random wave model in which 
the wave heights are described using the Rayleigh distribution. The obtained 
solution of the random wave field is used to predict the longshore currents. 
An empirical formula for determining the breaker parameters is developed, 
based on beach slope and incident wave steepness. The improved model is 
tested using an undulated bathymetry to validate the model physics. The 
model outputs of wave height and current are compared with data acquired 
from a wave tank and nat ieel beaches. The model is found to accurately 
forecast wave heights, breaker location, breaker type and longshore currents 
for several sets of conditions. Model limitations are discussed and 


recommendations for further improvement are made. 
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I. INTRODUCTION 


A. OPERATIONAL REQUIREMENT FOR WAVE FORECASTING 

Wave and surf forecasting has been an important part of Naval 
Amphibious operations since the beach landings of World War II. The need 
for the prediction of ocean and surf waves for military purposes led to the 
first -real attempt to quantitatively model ocean and beach waves. Ocean 
engineers and physical oceanographers, such as M.P. O'Brien, H.U. Sverdrup 
and W.H. Munk, made great strides in sea, swell and surf modeling during 
the 1940 - 1955 time period (Bascom, 1980). The wave and surf predictions 
made by the U.S. Navy are still based on these works. During the 
intervening time period, significant advances in knowledge of ocean wave 
dynamics have occurred, and it is time to update prediction techniques. 

Although warfare technology has advanced along all fronts, the standard 
amphibious beach assault is still a viable and preferable method of projecting 
forces ashore for many tactical scenarios. HQwever, amphibious assault craft 
are vulnerable to many types of wave dynamics. In general, landings are 
affected by the initial ocean sea and swell where the boats are launched as 
well as by the breaking waves in the surf zone, along with the accompanying 
longshore currents. Since in a major landing tens of thousands of lives are 
at risk, an accurate assessment of the sea, swell and surf is vital to a 
successful amphibious landing operation (Joint Surf Manual, 1976). 

In addition to the assault landing itself, an amphibious operation normally 
concludes with the construction of a temporary harbor facility’ to 


accommodate the resupply of the troops ashore. The survivability and design 


characteristics of these structures depend in large part on the environment 
around them. The ability to model wave height and direction, and longshore 
currents is necessary to operational planners and coastal engineers in 
designing and emplacing these facilities. 

In the past twenty years, separate models have been seoeigeed to deal 
with ocean waves, breaking surf and longshore currents. In 1981, the U.S. 
Navy contracted to have a model developed which would combine three 
complementary models. The goal of this comprehensive Sea, Swell and Surf 
Program (SSSP) model was to predict the wind generated wave height, 
calculate the surf zone characteristics based on the wind generated energy 
propagation, and use a two dimensional grid to compute the longshore current 
field. The model, developed by Wang and Chen (1983), was written in 
FORTRAN and designed to run on a main frame computer. To facilitate the 
model's use by Naval Oceanographers in the fleet, the model was converted 
into BASIC and implemented on the HP~-9845B/275 =micro-computer 
(Devendorf, 1985). The primary goal of this thesis is to modify the SSSP 
model to improve its computational efficiency. 

B. WAVE THEORY BACKGROUND 

The water surface elevation in the ocean can be viewed as a large 
number’ of individual wavelets at various frequencies. The waves that con- 
tribute the most energy to the ocean wave system are wind generated waves 
with periods of 1 to 30 seconds. The size and period of these waves are a 
function of the velocity of the wind, the duration of the storm and the 


distance over which the wind blows (fetch). 


The simplest theory is linear, small amplitude wave theory. Sea water 
is assumed to be incompressible and homogeneous, and surface tension forces 
are assumed to be negligible. The solutions of this theory describe the 
behavior of the individual components of ocean wind waves (Phillips, 1977). 
AS waves enter into shallow water, wave speed varies with the local water 
depth and refraction occurs. 

Initially, refraction of monochromatic waves was modeled using wave ray 
theory of a single wave (Arthur, 1949). It was not until the mid 1960's 
that the development of fast computers and the implementation of better 
physics enabled numerical analysts to study the generation of wind waves 
over an entire frequency spectrum. Early researchers in the wave spectral 
analysis field introduced dineetion variables to study straight and parallel 
contour wave spectrum transformation (Karlsson, 1969) and included the 
effects of bottom dissipation and local wind generation, 

As computers advanced, more complex refraction models were developed. 
Noda et al. (1974), using a relaxation finite differencing scheme, solved for 
the stationary wave spectral transformation. The model that is under study 
by this thesis was. originally developed by Wang and Yang (1981) and 
included the effects of bottom friction. Chen and Wang (1983) improved the 
model by including the effects of non-stationary waves. 

C. CURRENT MODELING TECHNIQUES 

Two numerical modeling techniques are currently used in wave and surf 
forecasting applications. The finite element method, while highly complex, 
has the advantage of being able to incorporate irregular localized bathymetry 
and has_ stable non-linear convergence properties. Wu and Liu (1985) 


introduced the method in their study of wave induced nearshore currents. 


10 


The finite difference method is a widely used numerical method. It is 

based on a linear differencing approach. Finite differencing normally uses a 
rectangular computational grid. This method was used by Noda et al. (1974) 
and Shaiu and Wang (1977) in their studies of nearshore circulation and wave 
energy transformation, respectively. The method is easier to implement than 
the finite element method due to its use of a regular grid describing the 
bathymetry, and is employed in the development of the SSSP (Wang and 
Chen, 1983). 
D. OBJECTIVES 

The objective of this thesis is to improve a sea, swell and surf program 
which is implemented on the HP 9845-B/275 micro-computer. The program is 
divided into subroutines, which include the open water wave generation 
module, wave transformation, the surf zone breaker calculations and the 
longshore current routines. To test the model, the output from the improved 
model is compared with data bases acquired from two natural beaches as well 
aS several wave tank experiments. A users manual has been prepared and an 
attempt has been made to interface graphic output with the current model. 
Limitations of the model are aiactiteea and future improvements suggested. 
This model is proposed for use by operational planners and Naval Ocean- 


ographers at sea. 
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II. OVERVIEW OF THE SSSP MODEL 


A. OVERALL PROGRAM DESIGN 

The Sea, Swell and Surf Program (SSSP) was developed by Wang and 
Chen (1983). The SSSP is essentially three complementary models which are 
merged to calculate offshore wave height fields, surf re information and 
longshore current velocities. The three modules , ocean, surf and longshore 
currents, can be run independently or the offshore wave energy can be used 
as input to the surf and current models. The SSSP is a numerical model 
which uses finite differencing to solve the governing equations. 

The original version, as written by Wang and Chen, was written in 
FORTRAN and designed to run on a mainframe computer. The version 
discussed in this thesis was converted to BASIC and implemented on an 16 
bit HP 9845-B/275 micro-computer (Devendorf, 1985). The micro-computer 
version consists of a control program with subroutines to calculate wave 
direction, height, and number, and output surf, swell and = current 
information. 

The SSSP considers an offshore wind generated energy spectrum which is 
used to estimate the wave growth based on initial inputs of wind speed and 
wave height or spectral energy bands. As the waves traverse down the user 
defined grid, changes in the deep water energy spectrum due to bottom 
friction and irregular bottom changes, are taken into account. If the surf 
zone model is to be run, a significant wave height, calculated from the area 
under the wave spectrum, and the peak frequency are used to provide a 


monochromatic wave height input into the surf model. It is also possible to 
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input monochromatic wave data directly into the surf module to compute the 
surf zone characteristics without running the computationally intensive 


offshore module (Wang and Chen, 1983). 


B. SSSP FLOWCHART 

The SSSP is controlled by one main program which, in turn, can call any 
of twelve subroutines. A general flowchart of the entire SSSP shows that 
the SSSP is divided into two main modules (Fig. 2.1). The user is prompted 
for all required information in a menu-driven, interactive fashion. Initially, 
the user chooses whether to run the offshore or the surf model. General and 
specific input is requested for the module being run. The open water or surf 
condition modules are calculated and the proper output is generated by either 
the OUTSW or OUTSRF routines. 

A more detailed flowchart of the open water condition module is 
presented as Figure 2.2. OPEN computes the wavenumber as a function of 
depth and then loops through a calculation series for each energy band. An 
adjustment is made for the angle of the swell/waves with respect to the grid 
axis. The DIRECT and HEIGH2 subroutines are called and the energy from 
each spectral band is summed and converted to wave height. 

The original SURFCON subroutine (Figure 2.3) used the most energetic 
energy band computed by the OPEN module as its frequency input, or used 
data entered directly by the user. The wavenumber was computed, the 
DIRECT and HEIGH2 subroutines were called, the average wave height was 
calculated and a series of breaker height and locations were predicted. 

If longshore currents were desired, the NEARCIR subroutine was invoked. 


The model, as implemented on the HP-9845, used an iterative process on a 
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Figure 2,1, General flowcnart of the Sea, Swell and Surf Program 
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Figure 2.2, Flowcharts of the Open water, wave Direction and wave 
Height subroutines (OPEN, DIRECT anc HEIGH2). 
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two dimensional grid to solve for the mean horizontal velocities. It was 
found that this scheme was computationally unstable as well as requiring as 
long as twenty-three hours to converge to a solution (Devendorf, 1985). 

To take advantage of the most recent advances in breaker height 
modeling, and increase the computational speed, the SURFCON and its 
ancillary subroutines were revised to use a one dimensional random wave 
model. A detailed flowchart of the new SURFCON module (Fig. 2.4), which 
replaces the original, shows the _ simplifications implemented by the new 
model. DIRECT2 is a simple Snell's law calculation, applicable over straight 
and parallel contours. HEIGH3 uses a random wave model which incorporates 
the probability density function for breaking waves. The output of the new 
SURFCON module includes breaker height, number of breakers, wave direction 
and surf zone width. 

LONSHOR is a new longshore current subroutine (Figure 2.4), which was 
implemented in place of the unstable NEARCIR to simplify and speed the 
calculation of the longshore currents. This model uses a simple radiation 
stress balance to calculate the current parallel to the beach. The next 
section will discuss the physics of the SSSP modules in detail. 

C. PROGRAM MODULES 
1. Offshore Module 
When the SSSP is run, a menu of initial choices is presented. If the 
mee chooses the "OPEN WATER" option, the OPEN subroutine is invoked. 
This offshore module computes the wave height and direction fields as a 
function of water depth, wind speed and direction, fetch, and the initial 
energy field. The initial energy can be described by a single significant 


wave height, direction and period, or as an implicitly entered energy 
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Figure 2.4, Flowcharts of NEWSURFCON, D4RECT2, HEIGH3 and LONSHOR 


subroutines used 


in the new model, 


spectrum. Water depth for the matrix can be read in from a data file or 
entered directly from the keyboard. A tidal height correction constant can 
be added to the depth field to take into account the state of the tide 
The bottom type can be specified, which determines the coefficient of friction 
for the wave height calculations. 
a. Wave Refraction 

Once al the required input is entered, the program begins its 

calculation cycle. The first field that is calculated is the wave direction 


matrix, using a steady state conservation of wavenumber (Phillips, 1977): 


3 3 : 
cae y= a = 7 
=x (k cos §) ay ‘k Sieriooy) 0 i) 
where: 
2 - 
k = wavenumber = = (m " 
L =' wavelength  (m) 
68 = wave direction (clockwise angle irom 
the X axis to the wave ray (see 
Figure 7) ) 
X =, offshore grid spacing (user specified 
(see Figure 5)) 
Y = longshore grid spacing (user specified 


(see Ficure 5)) 


The wavenumber is calculated, for the frequency band being considered, using 


the linear dispersion relation: 


Oo = (¢k tanh (k ny)i/? (2) 
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where: 


G = angular f¥equeney = 2ateecoen, 
Gg = acceleration of gravity (m sec *) 
h = water depth (1m) 

f£ = frequency = = (sec 7} 

T = wave period (sec) 


Since (2) is transcendental in k, a sixth order polynomial fit of (2) is used 
(Hunt, 1979) to provide a faster method of computing k than the iterative 
Newton's method used in the original Wang and Chen (1983) version of the 


SSSP. The wavenumber k is given as: 


k = ae + (] +0. 600675. - 0 : 
(gh) 


lv 
Ul 
S)) 
ui 
© 
—/ 


- 0.16084S7 “+ 0.063208 = 0.021748> 


+ 0.006548° + 0.001718’ + 0.000298° 


_, 1/2 
+ 0.000118”) +] 7 (3) 


where: 


tv 
Noakes 


Once the wavenumber is determined, subroutine DIRECT is 
called to compute the wave refraction of the wave components over the 
specified grid. DIRECT starts with initial values of the wave angles for each 


grid point, and uses a finite differencing scneme developed by Noda (1972) 
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to solve the differential equation for <4. .. The differential equation is 
] 


center differenced in the offshore (X) direction and forward differenced in 
the longshore (Y) direction. The finite difference scheme has greater weight 


on the forward differencing term to increase computational stability (Abbott, 


1979) and is given as: 


— : = ey 1 . a a \ 
aoe = sin K, os sin@s41/4- ae 2 ae) sins. 5? 





tat (is Sin @ 541 442) = sry (k cos ae 

= (k cos 8; 541) (4) 
where: 

ji .= Weichewms tactom = ga0e25 


A first guess for the wave angle is calculated using Snell's law 
solution for straight and parallel bottom contours. The numeric scheme is 
iterated until = Ht Cadet 0.005 of the previous calculation. The 10 X 13 
grid point system used by the SSSP is described by Figure 2.5. The initial 
differencing scheme used by the model is a simple forward stepping scheme 
(Fig. 2.6). The convention for the angle = is illustrated in Figure 2.7. 
Theta is measured with respect to a ray drawn perpendicular to the beach. 
Therefore, a normally incident wave would have a Theta value of 180°. 
After the direction matrix has been specified, the subroutine HEIGH2 is 


called to calculate the height transformation as the wave moves down the 


grid (towards the shoreline). 
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Figure 2.6. Numerical differencing scheme used by the OPEN water 
module, 


Shoreline 


a ee fT 


eo Ray 


iy 


Wave Crest 


Ray normal to beach 


Fiqure 2.7. Theta convention used by the model, 


2. Wave Height Transformation 
This subroutine takes into account the wave direction (previously 
calculated), the wind generation, and the bottom dissipation. HEIGHZ 
employs the Phillips-Miles growth mechanism for the wind generation 
calculations. An initial energy field is developed for each grid point by 


assuming that energy is conserved over straight and parallel contours: 


I 
> 
63 


sa Co cosmm 
Cox = 


where: 


tJ 
I 


wWaVe Cnergye (goules) 


group velocity (m sec 7) 


Cg 


The initial energy field is derived by integrating (5) to yield: 


G COSsa. 
So 


OCG Meese (©) 


where the subscript (0) is the initial condition at the offshore boundary. 
To simplify and increase the model's efficiency, Chen and Wang (1983) make 


four assumptions: 
a. Winds are steady 
b. Wave energy within” a “Spectral band 1s restricted tomimar 


band (Allows superposition of wavelets). 
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c. Each frequency component 


mean direction, 


om 


These four assumptions allow the steady state 


fi 
— 
af 


can be described by a single 


No wave-current interactions. 


energy flux for each spectral 


energy band to be written as (Longuet-Higgins and Stewart, 1960; 1962): 


0 


— a . = 
Paes oS recnGes =) + Wo —-(aeiryCo sin F) 
ON oy 
where: 
Sz) = “Spectral energy rlux 
oe Se 
o =—eGessiiy On Flurc {kom )} 
cg (=) = energy dissipation due to 
Beene 25 2LC LION 
E, (2) = energy generation cue to 


oD | 


The wind generation 


contributions: 


Wind stress 


= Meme anuwavwe Gl mee-10Nn 


term Ee , is the sum of the linear 


lo 
oo) 


and exponential 


(8) 


where: 


So = Phi Pa pee eve eee won mG. 
= ‘ 7 oy - = = 
So Illes Generation =a ocho. 


The energy generation term is based on the Phillips-Miles growth mechanisms. 
Starting with a flat ocean, wave growth is initially linear until larger waves 
are present, after which the growth is exponential. The linear growth term 
(Phillips, 1957) is based on random atmospheric pressure variations acting on 
a perturbed sea surface. The SSSP first calculates the pressure fluctuation 


term as a function of wavenumber and frequency: 


6 omer Vo a 
PAGKRe Ct) t= SS | SO 5 = 
_ ft _«= _ e _ a _ 
(9) 
where: 
W = wind speed in knots 
i) 
* 
eas 
vs = Ou. Se 
_ Gas 
e = Ores 
2 


>> 
iI 


angle petween winc and wave 


After P(k, Oo) is calculated, the linear growth term,a, is determined by: 


Z 3 
47 ke 
G. ad as aaa Pe Ga) (10) 


~ 


Pal 


where: 


. -3 
e = dénsity of water (kg m ~) 


The exponential growth term becomes important after waves have formed 
(LaBlond and Mysak, 1978). When flow separation occurs in the lee of a 
wave crest, momentum is transferred to the wave because of the pressure 


differential resulting in exponential wave growth (Miles, 1957; 1959a,b; 


1962): 
aoc W 
fo ee (=) cecma — (is) 
eel (a) b) 
where: 
aa —- Constant = pps20 
Spo “Sade O85 eircom wacer density 
=i 
C = wave celerity (m sec ~) 
Dee ecOnscanc. — 0, 90 
The Phillips-Miles wind wave generation has no mechanism for 
Shutting off the wave growth. In nature, wind waves will grow as a 


function of wind speed, duration and fetch. But after growing to a certain 
size and steepness, the waves will become unstable and break. Therefore, it 
is necessary to have a limiting function so that the waves will only grow to 
a "fully arisen" Sea. Several functions have been proposed to introduce an 
"equilibrium value”, above which the wave growth stops. The SSSP uses a 
Pierson and Moskowitz (1964) fully arisen sea in deep water. The high 


frequency tail is described by a Phillips (1957) equilibrium condition, 


bo 
=~) 


modified by Kitaigorodskii et al. (1975) and Thornton (1977) for intermediate 


or shallow water: 


SNe B gf "H(8)T(W.) (12) 


where: 





: 8 a, 

igh ue = 37COS 6 accounts for angular 
Spredcumd 
eS 2 jo (ee) ee 
a ae Co aaa = 5 ~ 
Sinh[2W 2 (W_)] 
n n 
W g (By 172 
n g 


2) meets the condition that ¢£ tanh(W I )=1 


Z 


T(W_) = 1 in deep water and for eee TC) —~ 1/2 Wand Siie = 


Gl 
2 -9 
1/2c0 = hedt z H( ¢ ) where q is on the order of 10 ~ (SSSP uses a value of 
0. 0720R 


Bottom dissipation, is calculated as the work done on the 


Eq? 
bottom due to the bed shear stress. Assuming the usual quadratic bed shear 


stress law, the wave energy dissipation for a particular spectral frequency 
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band is described by Hasselmann and Collins (1968): 


= 2 
€g(f) = oc, 1U,| ULC) (13) 
where: 
Cr = bed shear stress coefficient = 0.01 
U, | = total flow field bottom velocity 
U, (£2) = wave induced water particle bottom 
Velociey @rensd Specrrie Spectral 
component 


Using linear wave theory transfer functions to relate the surface elevation 
energy spectrum to the bottom velocity field, the frequency dependent 


dissipation function E468) can be expressed as: 
Se Ci aes 2) a0 | (14) 


= 7 he enone (kh)}. The total wave 


where the transfer function T(f) = [ g 
induced bottom velocity |, | is calculated by integrating the spectral 


transfer function across the entire frequency spectrum: 


oo iy 2 


Bottom dissipation in the SSSP is only calculated when deaing- with 
intermediate or shallow water waves. If deep water wave heights are being 


computed, the value ¢€, is set equal to zero (Wang and Chen, 1983). 
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3. Longshore Currents 

One of the most serious problems with the original model was the 
instability of the nearshore circulation calculations. As much as 23 hours of 
CPU time were required for the subroutine to converge to an answer. It 
was decided that a one dimensional longshore current model which assumes 
straight and parallel bottom contours would be substituted for the unstable 
two-dimensional model. This model uSes a simple alongshore radiation stress 
balance to calculate the mean longshore current. Although no circulation 
information is derived, the new longshore current subroutine, LONSHOR, 
provides the longshore current velocity efficiently and accurately. The 
Physics and details of the LONSHOR subroutine are outlined in Chapter Il. 

4. Surf Zone 

The original SURFCON module made calls to the same wave direction 
(DIRECT) and wave height (HEIGH2) subroutines as the offshore module. 
One disadvantage of this scheme is that these subroutines are _ very 
computationally intensive, especially when considering the narrow. width 
between grid points. A faster wave height model was needed. 

The state of the art in breaker height models uses a probability 
density function to describe breaking wave heights (Thornton and Guza, 
1983). Some limitations are imposed by using this new model, however. 
The model assumes that the beach has a local straight and parallel bottom 
contour. This bathymetric restriction allows longshore features such as 
sandbars to be modeled accurately, but offshore feature such as channels or 
canyons are smoothed out by a depth averaging subroutine. It is felt that 


this assumption is reasonable, in light of the fact that many beaches of 
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operational interest are of the "straight and parallel" type within the first 
two hundred meters from shore. 

The above assumption means that a simple Snell's law calculation of 
the wave direction is possible. A short subroutine, called DIRECT2, is 
implemented to calculate the wave direction matrix. The use of the new 
height and direction modules, HEIGH3 and DIRECT2 have increased the 
efficiency ‘of the SURFCON calculations immensely. The average run time for 
the surf zone model predicted wave heights has decreased from over four 


minutes to less than 30 seconds. 
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lil. MODEL IMPROVEMENTS 


A. OFFSHORE WAVE HEIGHT GROWTH 

The current version of the SSSP, converted to run in HP BASIC on the 
HP-9845B/275 micro-computer, gives reasonable values for wave height 
growth for wind speeds of 15 knots or less (Devendorf, 1985). However, for 
higher wind speeds of 30 knots, the SSSP exhibits a growth rate that is 
significantly faster than the growth of the JONSWAP spectrum (Bishop, 
1983), used for comparison in this analysis. The JONSWAP spectrum was 
chosen as a basis for comparison because the JONSWAP curves’ were 
experimentally derived from the limited fetch North Sea area. It is foreseen 
that the SSSP model may be applied to similar operational areas. 

For 30 knot winds, the model's wave growth not only occurs too 
quickly, reaching a maximum vane within 50 nautical miles, but the 
magnitude of the growth is too high (maximum wave height of 6.0 meters 
compared to a JONSWAP maximum of 4.2 meters. Before sensitivity tests 
were run on the wind generation portion of the model, the code was checked 
carefully to ensure that the equations were being implemented correctly. 

To determine the sensitivity of the model wave growth, several runs 
were made with changes to the key growth parameters. These parameters 
include A, the Phillips linear growth term, Equation 9 (Phillips, 1957); B, 
the Miles exponential growth term, Equation 10 (Miles, 1957; 1959a,b; 


1962), and R, a weighting factor set equal to 1.0 in the original model. 
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The first run was used as a standard to compare the modifications 
against. The inputs for the OPEN module were: 
Grid Spacing X = Y = 10 NM 
Wave Direction Theta = 180° (normal incidence) 
Wind (W) = 30 and 15 kts @ 270° 
oO 
Shoreline = 270 


Bandwidth (Df) 


0.01 Hz 


Frequency bands = 0.05 - 0.14 Hz 

The output for these "standard" model runs are shown as Figures 3.1 and 
3.2. The overbuilding of the open water waves for W = 30 Kts is evident 
when compared with the JONSWAP spectrum wave growth (Fig. 3.1). For a 
more moderate case where W = 15 kts, the model wave growth shows 
satisfactory agreement with the JONSWAP growth (Figure 3.2). To see the 
impact of changing the weighting factor, R, two runs were made setting R 
equal to 0.5 and 0.1 respectively. The predicted wave height with R equal 
to 1, 0.5 and 0.1. is shown in Figure 3.3. Setting R equal to 0.1 (chain 
dashed line) gives a ee fit for high wind speeds but underbuilds the wave 
field for lighter winds (Fig. 3.4). Clearly, a simple change of the 
weighting factor R will not resolve the problem. 

The Miles exponential growth term, B, is in the denominator of the 
energy calculation ty te value was increased to determine the growth 
retarding effect it would have. Two runs were made with B doubled (chain 
dotted line) and with B increased by a factor of 10 (not shown), depicted 
by Figure 3.5 The result of increasing the exponential growth term is to 
flatten out the upper portion of the growth curve which is not the correction 


that needs to be applied. 
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Because the model's rapid growth appears to occur in the linear region 
of the growth curve, the Phillips linear term, A, was analyzed closely to 
see if a simple constant multiple would bring the wave growth in line with 
the JONSWAP curves. A constant multiple of 0.125 is applied to the A term 
(Fig. 3.6). This method gives a better fit to the JONSWAP curve for high 
wind speeds, although the SSSP still over builds the wave heights. To test 
this multiple at lower wind speeds, a model run was made for a wind speed 
of 15 knots, for a bandwidth of 0.05 Hz and frequencies of 0.05, 0.10, 
0.15, and 0.20. The results of the model run are compared with the 
equivalent JONSWAP curve and the curve from the existing model (Fig. 3.7). 
The fit, while not exact, is reasonable, suggesting that the A term is 
responsible for the overbuilding of the model wave heights. 

The actual cause of the overbuilding of wind waves by the SSSP is 
unknown. Several ad hoc "fixes" to the model are suggested with the 
recommendation that further research on this problem be conducted. A 
simple constant multiple of the linear growth term, on the order of 0.10 to 
0.15 is one recommended empirically derived modification to the wind 
generation calculation. Another possible modification is to have the 
weighting function, R, decrease as the wind increases, to act as a damper 
on the over generation. At present, the SSSP wave generation routine 
remains in its original form. It is felt that the model is adequate for low 
wind speed generation. Operational planners should be cautioned about the 
high wind speed generation problems. 

B. WAVE REFRACTION MODULE 
To test the SSSP model's stability during the wave height prediction in 


the coastal zone, a sinusoidally varying bottom was explicitly entered using 
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the following equation; 


ane 
h = 0.025% + ch ()(1 - (4 ))cos(= =) (16) 
Cass b Ay He as 

where: 

E = bottom perturbation cCerm “(Cie 

: Ay = breaker depth (2.5 m) 
x = offshore distance to breaker line 
b 
(100 m) 
a = wave length of battom pertiumoamron (240 mm, 


This equation provides a smoothly varying bottom profile in both the xX 
(offshore) and Y (longshore) directions. To check that the bottom profile 
was being accepted correctly by the model, and that the calculated results 
were reasonable, approximately 25 model runs were completed with varying 
input parameters. Selected cases are presented herein for discussion. In 
the cases presented, the offshore waves are normally incident to the 
shoreline, Obliquely incident model runs were made with Theta ranging from 
+10 to -10 degrees from normal, for comparison with the normally incident 
case. 

1. Model Experiments 

The input parameters in the model are: significant wave height, H. 
(1 m); wave period, T (12 s); offshore boundary depth, h(12 m); artificial 
viscosity, T (0.0); and bottom perturbation term, £ (0.2). The beach 
faces 270° and the normally incident wave rays come from 270°. This gives 
a Theta value of 180° according to the model convention shown in Figure 


Loi 
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The model results for wave heights and directions are shown in 
Figures 3.873.195. The longshore depth profiles for locations inside (X=3) 
and outside (X=8) the surf zone are chosen to examine the wave refraction. 
Grid point X = 3 implies a distance of 40 meters offshore, well inside the 
breaker line, while X = 8 is 140 meters offshore and outside the breaker 
line. The depth profile presented in Figure 3.8 shows that the depth varies 
sinusoidally in the longshore direction. The beach profile along grid point Y 
= 3 is shown as Figure 3.9. As shown, the depths are slightly deeper than 
the plane beach inside 100 meters (breaker position), and shallower outside 
MOG Bmeters , 

AS waves enter the surf zone, their increase in height is a function 
of water depth. For the points outside the breakers, the largest waves are 
expected to occur where the water depth is the shallowest due to shoaling. 
As can be seen from Figure 3.10, the peak waves (dashed) are towards the 
edges of the domain, but they do not coincide with the minimum in depths at 
0 and 240 meters. This anomaly is a result of the boundary conditions of 
the model, Also note the asymmetric peaks at 40 and 200 meters. This 
implies errors in wave height prediction due, in part, to the wave direction 
calculations. These errors will be addressed later in this discussion. 

Wave direction is calculated in terms of Theta. An angle greater 
than 180 tends to turn the wave ray to the left while Theta less than 180 
turns the ray to the right (See arrows in Figure 3.11). The model at least 
quantitatively predicts the wave refraction pattern correctly for all cases, 
consistently turning the rays to “attack" the shallower water depths outside 


the breaker line. The model's boundary conditions force €¢ fe ¢. and 6 


Z ne 1 


= § - where 6 is the wave direction at the a long shore grid point. This 
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accounts for the departures of the wave direction at the edges of the 
domain. The dotted lines are the interpretation of how the wave direction 
patterns should behave. In the middle of the domain, where both cases 
exhibit a Fora point, the wave direction passes through 180 or normal to 
the beach as expected. 

The final test conducted compared the effect of Tau on the wave 
direction calculations. Tau is an artificial viscosity parameter which is 
intended to be a smoothing term for the calculation of wave angle. 
Changing Tau from 0.00 to 0.25 had a negligible effect on the wave height 
calculations, but the difference in direction increased with decreasing 
distance from shore. The case where X = 8 is plotted in Figure 3.12. The 
introduction of the Tau term (dashed) decreases the magnitude of the wave 
angle, but the maximum (minimum) values are not exactly in phase with the 
Tau = 0.00 case. The maximums are shifted slightly toward the edges of the 
domain. This is the expected behavior of the smoothing term. 

The dissipative model exhibits unexplained behavior inside the 
breaker line as depicted by Figures 3.13 - 3.15. At X = 5 (Fig. 3.13), 
the Tau model reduces the magnitude of the wave direction field. At X = 4 
(Fig. 3.14), spurious asymmetry appears in the direction calculations with 
the Tau term applied. Furthermore, at X = 3 (Fig. 3.15), the effect of Tau 
is to increase the magnitude of the wave direction. This apparent angular 
deviation has no effect on the model's wave height calculations because, 
inside the breaker line, wave height becomes a function of depth only. 
However, this instability does effect the longshore current calculations, 


which are sensitive to local wave direction inside the surf zone, 
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Another indication of a numeric error in the wave direction 
calculation was noted in the exact center (Y = 7) of the domain. Since the 
depth profiles exhibit a maximum (minimum) value at the domain center, 
theory Brecests that the waves remain normally incident at these nodal 
points. Therefore, the Theta value should be 180 degrees at the center of 
the domain. The model output, both with and without Tau, exhibits a 
monotonically increasing error as seen in Figure 3.16. 

2. Improvement of the Numerical Scheme 

In an effort to eliminate the numerical instability in the wave 
direction module, the finite differencing scheme was analyzed to see if an 
improvement could be made. The original scheme, pictured in Figure 2.6, 
uses a simple centered and forward differencing method and is fairly robust in 
the longshore (Y) direction but is less so in the offshore (X) direction. To 
increase the accuracy of the numerical scheme, the method pictured in 
Figure 3.17 was. substituted in the DIRECT subroutine for the _ original 
scheme, This particular method is a well accepted technique which has 
proven to be stable and free of damping (Haltiner and Williams, 1980). 

Several model runs were made to test the new technique. Using the 
standard input parameters and setting Tau equal to 0, the first result was a 
gain in accuracy of the wave direction calculations. The domain's central 
directions were exactly 180° (to three significant figures), as seen in Figure 
3.18, and the refraction angles were symmetric about the center of the 
domain as expected due to the symmetrically varying bathymetry (Figure 
Opel ol). 

A second result of the implementation of the improved numerical 


technique is the stabilizing effect on the wave height calculations. The 
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previous asymmetric wave heights are now symmetric with respect to height 
and peak location (Fig. 3.20). 
3. Boundary Conditions 
The SSSP model imposes no flux boundary conditions on_ the 
refraction and wave height calculations so that the left and right boundaries 
are specified as follows: 
Q(M,1) = Q(M,2) 
Q(M,N) = Q(M,N-1) 
where: Q = the quantity being calculated 
M = offshore row number (10 to 1) 


N = the last longshore column (13) 


The offshore row is completely specified by the initial conditions. The 
calculations then proceed from offshore to onshore and from left to right on 
the grid. 

No flux boundary conditions account for the unusual behavior of the 
wave height and direction curves at the edges of the domain. To satisfy the 
condition of well-posedness for the applied model, cyclic boundary conditions 
are substituted such that; 

Q(M,1) = Q(M,N-1) (17) 
Q(M,N) = Q(M,2) 
A plot of wave height versus longshore distance at X = 8 using the same 
inputs as the previous case is presented in Figure 3.21. Note the shift of 
the peaks toward the edges of the domain as well as the abrupt increase in 
the magnitude of the wave height. These numerical results now match quite 


well with theory which requires that the wave rays turn toward the 
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shallowest water, implying that the largest wave heights should occur at the 
edges of the domain, 

The wave direction calculations also show improvement. Close 
agreement between the cyclic boundary condition values and the expected 
values (dotted) is shown in Figure 3.22. Additionally, the average value on 
the boundary segment agrees closely with the exact value (solid circles). 

4. Results 

Based on the results of many model experiments, the new model 
appears to offer an increase in stability and computational accuracy. The 
improved wave direction subroutine consistently turns the wave rays in the 
correct direction. The magnitude of the direction changes is qualitatively 
correct. The wave height subroutine works well both inside and outside the 
surf zone. Wave heights are highest at the shoal and lowest at the trench. 

One drawback to the original wave height subroutine is the 
computational speed. When running the program to generate an energy 
spectrum for sea and swell, the model takes approximately 45 minutes per 
frequency band. Even the monochromatic surf module is slow (4.5 - 5 
minutes) because of the number of iterations required to converge to the 
preassigned criterion. 

Least acceptable of all is the Nearshore Circulation module which 
proved to be unstable as well as requiring as much as 23 hours to 
asymptotically approach an answer. To provide fast, accurate longshore 
current information, a one dimensional model is substituted for the two 
dimensional circulation model. While some potential rip current information 


will not be provided by this simple model, the one dimensional model has the 
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advantage of being fast, accurate and able to provide the most important 
input, the longshore current distribution. 

A state of the art, one-dimensional wave height model is used as 
input to the longshore current model. The new wave height and longshore 
current subroutines are described below. 

C. ONE DIMENSIONAL WAVE HEIGHT CALCULATION 

To increase the SSSP's computational efficiency while maintaining 
reasonable accuracy, a one dimensional random wave height model was 
substituted for the original two dimensional monochromatic model. This 
model characterizes the transformation of the wave height probability density 
function (pdf) from offshore to the shoreline (Thornton and Guza, 1983). 
The model assumes that the waves are narrow banded in frequency and 
direction such that the waves are specified by a single mean frequency f and 
mean direction 6, that the bottom contours are straight and parallel, and 
that the wave conditions are stationary. This probabilistic approach to the 
wave breaking problem was originally suggested by Collins (1970) and Battjes 
(1972). The form that is used here is based on integrating the energy flux 
balance equation. The calculated wave heights are a function of the shoaling 
conditions of the integral path from offshore to the shoreline (Battjes and 
Janssen, 1978). 

For straight and parallel contours, the energy flux is balanced by wave 


breaking and the frictional dissipation: 


(CE CGmeOsuc) (18) 
a = oe” + _oe- 
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where: 


> 


Se) dissipation due to wave breaking 


frictloneameaissipazmen 


<p> 
In the domain of interest, it has been shown that the frictional dissipation is 
less than 3% of the dissipation due to breaking waves (Thornton and Guza, 
1983) and will be neglected. 


The waves in the new model are described by the Rayleigh wave height 








distribution: 
2 
a 
H (19) 
2H rms 
p(H) = 7 = 
rms 
where: 
H = root mean square wave height 


INS 


where aS is the root mean square wave height. This distribution takes 
into account the fact that there is a distribution of wave heights throughout 
the entire field for a given frequency. Previous random wave models 
truncated the Rayleigh distribution in various ways. Recent observations 
have found that even after the waves break, the Rayleigh distribution -is still 
valid (Thornton and Guza, 1983). Therefore, in the improved model, the 
waves are everywhere Rayleigh distributed. The probability of a wave 
breaking Pp, (H) is described by a _ simple weighting of the Rayleigh 
distribution: 


p, (H) = W(H) p(H) 
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The weighting function, W(H), which makes larger waves at a given depth 


more likely to break, is given by (Thornton and Guza, 1983): 


H 
W(H = iS 
(H) es, (20) 
where: 
Y = adjustable parameter based on 


field data 


The model assumes that breaking waves are similar to periodic bores 


(LeMehaute', 1962) as depicted in Figure 3.23. The bore dissipation per 


unit area is calculated: 


E Leas ane z lo, (BH) > 0 = 
y « “SReGe foe wai —ooT an 
bore 4 n,n, 4 ne 
where: 
QO = volume discharge across bore 
B = breaker coefficient (described 


from Gata) 


The volume discharge parameter is described for a linear periodic bore 


(Hwang and Divoky, 1970) by: 
Q = Ch/L (22) 


where C is the wave celerity. Using linear wave theory, the wave energy 
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which describes how much of the wave is breaking at any 


Periodic bore used to describe spilling breakers (after 


Thornton and Guza, 
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time, 


given 


and group velocities are given by: 


Ee sae H H) qQH = =, 
g J p (H) g 09 He (23) 
and, 
6 cee 2kh 
‘ 8x Di actaiae ie 


(24) 


After integration and the substitution of the probability density function, 
Equation 21 leads to: 


a 


=e 
ae re 
= gt og -p 





where: 


(25) 


f = 


mean frequency 


This equation for the bore dissipation can be tuned to the data by varying 


the breaker coefficient B, and the adjustable parameter 


Y . Substituting the 
bore dissipation function (25), the energy flux equation (23) is given by: 

ul 2 

= C 

(5 Og a 


nae Coame.) 


sa 
wa 


. ee £B 7 (26) 
me su ee ee ae eae 
ae 
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For shallow water, the group velocity Cg is approximately equal to the wave 
celerity C which equals (gh) Equation 
the wave height module because, 


(26) was chosen to be used in 


although a more complicated expression 
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gave a Slightly better fit to the data, Equation (26) provides an analytical 
solution which allows a detailed study of the model's’ behavior. The 
analytical solution is obtained for normally incident waves impinging on a 
plane sloping beach and is used to compare the accuracy of the numerical 
model. The analytical solution is solved by integrating (26) for the wave 


height (Thornton and Guza, 1983): 





Ase) 10) ee sy 3. _ a =175 
ee ah [Ao aaa "2374 ee (27) 
O O 
where: 
4 
_ 2390.72 \ eaters 

; Et oS a 
Bae 

tan 8 = slope of beach 

y =f pee 

O emetic. 


The numerical model uses a modified Euler integration scheme to 
integrate the energy flux equation. A weighted predictor/ corrector loop is 
used to increase accuracy and stability as the water depth becomes shallow. 


The basic calculation is a forward stepping scheme with a variable grid size: 


Mt 
B Co = ecg + =(<e, > + <e, >} ) AX 
Well outside the surf zone the wave height changes are very slight with 
distance as the bore dissipation is essentially zero, and are relatively 
insensitive to the grid size. However, once the waves start to break, the 


bore dissipation increases rapidly, modifying the wave heights. To increase 
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the resolution, the mesh size decreases by a factor of four when the bore 
dissipation reaches a value of 0.01. This flag value was arrived at by a 
series of model tests under different initial conditions. 

The model was tested against the analytical solution for a variety of 
initial values. The initial conditions insured that shallow water criteria were 


1/2 rather than 


met so that the shallow water approximation of Cg = (gh) 
the explicit Equation (24), could be used. The first test used a shallow 
steepness offshore wave (Ho /Lo = 0.01) having a height of 1.0 meter 
normally impinging on a plane beach with a slope of 0.04. The period of the 
wave was 10. seconds. The numerical and the analytical solutions are 
compared in Figure 3.24. Agreement is very close through most of the 
offshore grid points. As the water depth gets very shallow, the step size is 
again reduced to insure the solution converges for the last few grid points. 

For a steep wave case (H/L, = 0.02), an incident wave of 2.5 meters 
was used as the offshore input (Fig. 3.24). The wave period, beach slope 
and normal incidence all remained the same. The numeric and analytic 
solutions are again in good agreement, although the numeric solution slightly 
underpredicts the wave height as it did in the case where Ho equals 1.0 m. 
It is of interest to note that the high wave steepness case shows no growth 
in wave height due to shoaling prior to breaking. The wave heights decay 
monotonically toward the shoreline. 

Having satisfactorily compared the model results with the analytic 
solution, further improvements to the model's computational efficiency were 
made, Since the model is one dimensional, the depth variations in the 
longshore direction are averaged out. Sand bars and other longshore features 


are taken into account but offshore trenches and channels are smoothed out 
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by this process. The new height module, with this simplification, runs six 
times faster than the original module. 

The new wave height subroutine has many advantages. It uses a simple 
but robust probabilistic approach to the wave height problem. It allows 
variable bathymetry in the offshore direction and obliquely impinging wave 
fronts. It is fast, accurate, and the governing equations can easily be fine 
tuned by’ varying the parameters. 

D. LONGSHORE CURRENT CALCULATION 

The original longshore current model (Chen and Wang, 1983) used a two 
dimensional momentum flux balance which was numerically intensive. The 
model converged very slowly (over 23 hours of CPU time in some cases), and 
the calculations exhibited unstable behavior (Devendorf, 1985). To provide 
the necessary longshore current information to an operational planner, a 
simple one dimensional model was substituted for the original two dimensional 
subroutine. The model was designed to work with the one dimensional wave 
height model discussed in the previous section. The random wave probability 
model is used as input for the longshore current calculations. A narrow 
banded wave distribution is assumed so that the wave frequencies can be 
approximated by a mean frequency f and direction = . The longshore 
momentum equation is a balance between the mean wave momentum flux 
divergence and the longshore component of the bed shear stress (Thornton 


and Guza, 1985): 


(29) 
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where: 


S an time and depth averaged covariance 
Y between unsteady velocity components 
b 
ty = bed shear stress in longshore direction 
s is made up of two terms which are assumed to be _ statistically 


yo. 


independent. The first term is the radiation stress term (S.) which is the 
wave induced momentum flux. The second term (Soy is the depth 
integrated turbulent Reynold's stress, often parameterized in terms of an 
eddy viscosity term. 

Monochromatic current models require the eddy viscosity term to smooth 
out the discontinuity in the current field at the breaker line. Since no 
waves break outside the surf zone in a monochromatic model, there is no 
current generated. However, there is a jump discontinuity at the breaker 
line which requires smoothing by the eddy viscosity term. In the random 
wave model, there is no clearly defined breaker line. Some waves break 
offshore and, as the beach is approached, more waves break until, in the 
inner surf zone, almost all waves are breaking (Thornton and Guza, 1985). 
The smooth transition of energy dissipation reduces the requirement for the 
eddy viscosity term, which simplifies the velocity calculation. 

The radiation shear stress Sy is described by (Longuet-Higgins and 
Stewart, 1964): 
=) S ee 


yx 7 EC. cos 6 


0 
C) 


(30) 
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Snell's law requires that: 


_ sin 6 
ee one fant oe 
C Co 31) 


so that when (30) is differentiated the result is: 


si Sane Sin 86 


= O <4 
ax, Sox = acme -(E ae cos 9) = [el cao 
O 


er 





= 


ae 
ms 


(32) 


This implies that the shear stress is balanced by the bore dissipation for 
breaking waves on a straight and parallel beach. 


The bed shear stress is given as: 


i, eee] (33) 


where: 


instantaneous total bottom velocity 


cy 
i 


< 
tl 


longshore current velocity 


The bed stress is linearized by assuming small angles of wave incidence and 


weak longshore currents (Thornton, 1970; Longuet-Higgins, 1970): 


ie m= VSS = |b Nt (34) 


Using shallow water wave theory, the bottom velocity is computed by: 


oo 1/2 - 2 
<i ee! vt 2 
HE ce Saye) | ac auaie le ale 
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Making the appropriate substitutions into the longshore momentum relation and 


solving for V: 





6 
Oa. ge B> Fg’? sin 6 Mrms 
: Cey" Cy 7 (36) 


The numerical scheme is compared to the analytical solution, which is 
obtained by simply substituting the expression for Hea (Equation 27) into 
the velocity calculation (Equation 36). The numerical and the analytical 
solutions for the shallow and steep wave cases studied above are compared in 
Figures 3.20 and 3.26. The numerical solution of V_ slightly underpredicts 
the current velocity compared with the analytic solution. This is due in 
part to the introduction of a small error in Ba which was discussed in the 
previous section. The value of os is raised to the sixth power in Equation 
(36) which means that V may be sensitive to small errors in H aes 
However, the magnitude of the error for V is on the order of 0.05 and 0.07 
m/sec which may be neglected. 

The initial runs of the updated wave height and longshore current 
subroutines are presented and the reworked SURFCON module is discussed in 
the next section. 

E. IMPROVEMENTS TO THE SURFCON MODULE 

During the reprogramming effort, a major consideration was maintaining 
the integrity of the original model, while incorporating modular 
improvements. The new height, direction and current subroutines were coded 
under new names, leaving the old subroutines intact. To maintain 
consistency, the control module, SURFCON, was renamed OLDSURFCON and 


a new module, NEWSURFCON, was written to make use of the new 
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subroutines, The advantage of this approach is that the model can be 
changed from the new to the old version simply by changing the main 
program "CALL NEWSURFCON" to "CALL OLDSURFCON". 

The NEWSURFCON module is basically the same as the old module with 
some exceptions. The wavenumber k, celerity C and group speed Cg are 
initially calculated. Then the new wave direction and height subroutines are 
called, which return direction, height and longshore current matrices. 

The first major change incorporated in this new module is the way the 
breaker line is defined. In the old module, a quadratic equation was solved 
for the location of the first breaker line. Since the new breaker heights are 
calculated using a weighted probability density function, Pp, (H) = W(H) p(H), 
there is, by definition, no specific spot where all waves are said to be 
breaking. The weighting function, W(H), as defined in Equation (37), gives 


a measure of the percent of breaking waves (Thornton and Guza, 1983): 





H 4 
WH) = [==] =e Oe (HH) (37) 
gla 0 b 


To specify a "breaker line", the model first differentiates between steep 
waves, which show no wave growth due to shoaling towards the shore, and 
shallow waves, which increase in wave height until they break. The critical 
wave steepness, determined from Equation (27), is Ho/L = 0.02. A wave 
with a steepness greater than this criteria is considered steep while a wave 
with steepness less than 0.02 is considered shallow. For steep waves, a 
decision must be made as to what percentage of waves having broken 
constitutes the surf line. After some rather ad hoc experimentation, a value 


of 0.33 was chosen as a reasonable, which in effect defines a "significant 


breaker height". That is, once 33 percent of the waves are breaking, the 
model chooses that point along the grid as the first breaker line. This test 
is actually conducted in the wave height subroutine. As the model moves 
down the column, the parameter W(H) is tested at each grid point to see if 
it greater than or equal to 0.33. Once this value is reached, the breaker 
height, breaker depth, and distance from shore are returned to _ the 
NEWSURFPCON module. For shallow waves, the model selects the maximum 
wave height as the breaker line. These two methods give good estimates of 
the observed breaker line, as discussed in Chapter 5. 

Once the initial breaker height is defined in this manner, the number of 
additional lines are computed by dividing the initial breaker distance by the 
wave length corresponding to the average surf zone depth. If there are one 
or more wave lengths between the initial breaker line and the beach, the 
number of breaker lines are incremented accordingly. No additional 
calculations are performed on these interior breaker lines, if present. The 
breaker type is calculated based on the value of the surf parameter, Eta, 
which is a function of the beach slope and the deep water wave steepness 


(Battjes, 1974); 


E tan & a tan 8 
eZ H (38) 
(HQ/L) (Qn ( Q),1/2 
Cet 


If Eta is less than or equal to 0.4, the breaker type is defined as spilling. 
Eta values between 0.4 and 2.0 are defined as plunging and values greater 


than 2.0 are called collapsing or surging waves. 
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The next calculation performed by NEWSURFCON is the "Effective surf" 
computation. Effective surf is a measure of surf intensity and is used by 
operational planners to provide a criterion for the feasibility of conducting 
amphibious ‘landing operations using specific kinds of equipment. Effective 
surf, expressed as a wave height in feet, is a modification of significant 
breaker height, taking into account the total beach and surf conditions. It is 
important to note that effective surf is simply a planning parameter, and not 
a direct correspondence to the significant breaker height. Most commonly, 
H ee is larger than H. but in some cases, especially when the H. is small 


and the breaking waves are all spilling, H can be smaller than H.- 


eff 

The calculation of effective surf uses a simple "look-up table" method 
which is outlined in the Joint Surf Manual (COMNAVSURFPAC/ 
COMNAVSURFLANT Inst. 3840.1, 1976). The calculations of effective surf 
are retained from the original model with only minor’. modifications. 
Originally, the effective surf was retained as a wave height matrix. 
However, H ee is simply a single number for a specific beach, based on 
significant wave height and direction, period, wind speed and direction, and 
percentage of spilling waves. The new module calculates a single number for 


H and outputs all the Hy parameters as specified by the Joint Surf 


eff ff 


Manual, Other information of use to the planner is printed at the same 
time. The new calculation also allows for a wind input (speed and 
direction) which was not available previously. 

The final operation in NEWSURFCON is the limitation of the wave heights 


to the maximum stable wave height for a given depth and period. This 


a 


maximum wave height is given as (Wang and Chen, 1983): 


zal 


H = 0.067 tanh (kh) ( —) (39) 


max 
Each generated wave height is compared against this maximum value and is 
reset to Bae if it exceeds the limiting wave height. The wave height, 
direction, current and effective surf matrices are then sent to OUTSURF for 
formatted output to either the screen or to a designated printer, 
F. DEVELOPMENT OF AN OPERATIONAL USERS GUIDE 
A simple users guide was written for the new version of the SSSP (Gill, 
1985). Several guidelines were adhered to during its development. The 
manual assumes some familiarity with the HP-9845/275 and a_ working 
knowledge of oceanographic terms. The manual is tutorial in nature. It is 
divided into two parts, for the two main modules OPEN and SURFCON. The 
manual explains each input request and appropriate examples are included 
where necessary. An effort was made to keep the manual brief and user 
friendly. The model, being extensively menu driven, is straight-forward 
enough to require little documentation, after it has been used once or twice. 
The BASIC code is extensively documented throughout. Where appropriate, 
the old program lines are left in place (commented to prevent execution), 
with the replacement lines documented to point out the reason for the 


specific change. 
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IV. DISCUSSION 


A. SENSITIVITY TEST OF BREAKER PARAMETERS 

As discussed in Chapter III, the computation of : using the Rayleigh 
probability density function model, depends on two breaker parameters, 
Gamma and B. Gamma is an adjustable coefficient of O(1) which appears to 
be strongly related to the beach slope (Sallenger and Holman, 1985). 8B is 
the breaker coefficient which represents the percentage of foam on the face 
of the breaking wave (i.e., the intensity of the breaking wave). This 
parameter is expected to be less than or equal to 1.0, since a B value of 
1.0 corresponds to a fully developed bore (a wave that has fully broken from 
crest to trough). Values of B greater than 1.0 imply that the bore 
dissipation function underestimates the wave dissipation due to breaking 
(Thornton and Guza, 1983). 

To verify the validity of the model assumptions, 24 wave data sets were 
studied to gain an understanding of the interaction of the B and Gamma 
values as used in the model. These data sets are summarized in Table 4.1. 
The data were digitized and appropriate Gamma values chosen, based on the 
slope of H ams! Dt inside the surf zone. The Gamma values ranged from 0.39 
to 0.60 and, for each data set, the specific Gamma value was held constant 
while the B values were varied. A curve fitting routine was used to 
optimize the B parameter for a given Gamma, As B was varied, the new 
model curve was compared with the real data points, and a sum of the 
square of the errors was generated. The B values were iterated until the 
error sum was minimized. Plots of the model output versus the data points 
were generated for each case, 
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A total of 24 data sets were tested against the model (See Table 4.1). 
In every case, the model accurately predicts the wave heights described by 
the data points. The first ten data sets are from a laboratory study of 
energy saturation of irregular waves during shoaling (Vincent, in press). 
These wave tank tests were conducted with a plane bottom slope of 1:30, 
with various frequency and incident wave heights. The experiments can be 
divided up into two groups based on wave steepness. The low steepness 
cases OM < 0.02) are represented by Case 1638 (Fig. 4.1). The 
Gamma value of 0.60 was derived from the slope of H ems /h as described 
above. Using the curve fitting routine, an optimum B value of 0.82 was 
chosen. Note the good fit of the model curve (solid) compared with the 
data points (circles). The high steepness cases are represented by Case 
1148 (Fig. 4.2). Gamma was set equal to 0.60 and a B value of 1.2 was 
obtained. As before, there is good agreement between the model curve and 
the data points. 

To determine the optimized model's behavior when used with more 
difficult bathymetry, data from several wave tank experiments with offshore 
bars were tested (Battjes and Stive, in press). Case 5 (Fig. 4.3) has a 
simple offshore bar with a mean bottom slope of 1:25. Gamma was estimated 
to be .45 and a B value of 1.2 was obtained. Although the presence of the 
offshore bar (Fig. 4.4) radically changes the wave height data profile, the 
optimized model curve again faithfully predicts the wave heights. Wave tank 
test case 15 (Fig. 4.5) has a very complex bathymetry (Fig. 4.6). The 
model divides the on/offshore distance into 10 grid points with the 


intermediate depths calculated by straight linear interpolation. Therefore, 
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complex bathymetry, such as that in Case 15, was smoothed out to a 
certain degree, Even with that limitation in mind, the optimized model 
curve appears to accurately estimate the recorded wave heights. The model 
results will’ be tested against field data in the next chapter. 

It is evident that the random wave model, with the appropriate choice 
of the surf parameters, Gamma and B, will accurately forecast wave heights 
throughout the domain. A _ sensitivity test for the choice of parameters was 
conducted on three representative cases 1638, 1148, and 5. Case 5 (Battjes 
and Stive, in press) was chosen for discussion here because of its more 
interesting bottom profile. It was felt that this barred profile was 
representative of many real beaches and would provide a more rigorous test 
than a simple plane bottom. Using a Gamma value of 0.45, the optimum B 
value of 1.3 was changed to 1.7, 1.5, 1.1 and 0.9 (Fig. 4.7). B values 
less than optimal tend to over-estimate the wave heights and vice-versa. It 
is of interest to note that changes in B of 22-33 percent only results in an 
increase in the model error of about 10-12%. This was confirmed in a study 
of Torrey Pines, California wave data by Thornton and Guza (1983). They 
found that a variation of +/- 25% about the optimum B value resulted in an 
increased model error of less than 10% (Fig. 4.8). 

It was found that the model demonstrated a similar sensitivity to changes 
in the Gamma value. A change of about +/- 20% to the optimum Gamma 
value of 0.45 resulted in an increased model error of 8-10% (Fig. 4.9). 
Sensitivity tests were conducted on high and low steepness plane beach cases 
(1638 and 1148) with similar results. In the next section, techniques and 
rationale for choosing B and Gamma values, based on initial conditions, are 


discussed. It is important to bear in mind the above sensitivity analysis. 
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The model appears to be fairly robust with respect to small changes in the 
surf parameters. Therefore, slight errors in the choice of these values 
should not degrade the final performance of the model to a great extent. 

B. CHOOSING THE BREAKER PARAMETERS 

The two breaker parameters, B and Gamma, are multiplicative constants 
in the Hae wave height calculation (see Equation 27). They could 
therefore, be combined as a single bulk modulus coefficient, or kept as two 
separate parameters. It was felt that keeping the two terms separate would 
better serve the requirements of the model. 

Some work has already been done on the problem of choosing the breaker 
parameters. A recent study by Sallenger and Holman (1985) has determined 
a functional relationship between Gamma and the average beach slope based 
solely on field observations inside the surf zone. This relationship is given 
by: 


y = eee ange 
(40) 


Using the 24 wave data sets, the chosen Gamma values are plotted against 
beach slope (Fig. 4.10). The solid line is a straight linear regression of the 
data points while the dashed line is the Holman/Sallenger function. There 
seems to be a general agreement in the slope of the line although tne 
intercept is slightly different. Gdmma is also plotted against wave steepness 
(Fig. 4.11) but the plot is less linearly correlated. As an initial trial, the 
Holman/Sallenger functional relationship will be adopted as the model's 


convention for choosing the Gamma value. The average slope is defined here 
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as the depth at the outermost grid point divided by the offshore grid 
distance. 

Choosing a proper predictor for the B value proved to be more difficult. 
Since B represents a measure of breaker intensity, it is felt that it would be 
related to incident wave steepness which has previously been shown to be 
related to breaker type (Battjes, 1972). The optimum B values for the 24 
wave sets are plotted against steepness in Figure 4.12 and against beach 
slope in Figure 4.13. The wave steepness in this case was defined as the 
incident wave height divided by the deep water wave length (Ho /L,). 
Although there is a considerable scatter, B appears correlated with wave 
steepness and less with beach slope. 

Several other parameters were considered as predictors for B. The B 
and Gamma values were compared with the surf parameter, Eta = tan 8 / 
(ge alae Additionally, the Bulk Modulus (B°/Gamma‘) was plotted 
against beach slope and wave steepness. Since, from Figures 4.11 and 4.12, 
B appears to be related to both beach slope and wave steepness, it was 
hoped that Eta, which is a measure of the expected breaker intensity, would 
offer a predictor for B. The plot of B versus Eta (Fig. 4.14) shows an 
interesting pattern. B appears to increase almost linearly as Eta increases 
to about 0.6 and then decreases in an exponential fashion to an asymptotic B 
value of 0.8. This behavior appears reasonable. As the surf parameter 
increases, indicating a transition from spilling to plunging waves, the B value 
(amount of wave face that is breaking), should increase. 

The maximum value that B should reach, according to theory, is 1.0. 
Studies have shown, however, that in many cases, the optimum B value is 


on the order of 1.4 (Thornton and Guza, 1983). The maximum B value 
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occurs when Eta is about 0.6, indicating that the breakers are mixed spilling 
and plunging. As Eta continues to increase, the B values fall more into line 
with theory and approach the expected value of 1.0, as the waves become 
fully plunging. 

_ To simplify and expedite the model development, the following piecewise 


continuous function is used to determine the B term (Fig. 4.14): 


¥ 


B= 2.167 Eta + 0.4 Eta <= 0.6 (41) 
B =- 1.70 Eta + 2.7 0.6 < Eta <= 1.0 
B= 1.0 Eta > 1.0 


This choice is somewhat ad hoc and is based primarily on the data sets 
discussed above. Future refinements in the techniques for predicting the 
breaker parameters may. substantially improve the use of this model 
operationally. Using equations (40) and (41) as predictors for Gamma and 
B, the SURFCON module is essentially complete. The next chapter will 
discuss testing the model against two well-documented natural beach data 
sets, to get a feeling for how the improved model behaves against actual 
data, 
C. MOTIVATION FOR THE RANDOM WAVE MODEL SUBSTITUTION 

As discussed in Chapters 3 and 4, the primary motivation for substituting 
the one dimensional random wave model for the two dimensional model was to 
increase computational efficiency and accuracy. The improved model 
calculates longshore current information quickly and accurately. Significant 
wave height predictions compare well with actual data sets, as will be 


discussed in Chapter 5. There is , however, a potential for loss of detailed 
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information by using the one dimensional model. Nearshore circulation and 
longshore variation in wave heights are smoothed out by the one dimensional 
approach. 

The surf model (both old and new versions) uses a single wave height 
input to define the offshore boundary condition. This implies that variations 
in the interior wave heights and directions are due only to variations in 
bathymetry. Because the model uses a very small grid (only 13 longshore 
grid points), very accurate bathymetry must be entered. It is questionable 
whether the operational planner would have access to bathymetric data of 
the required accuracy to define a two dimensional bathymetry. 

Many beaches which would potentially be used for amphibious operations 
are shallow, almost planar in nature. Other common landing beaches are 
barred, but rather uniform in the longshore direction. Therefore, for many 
operationally relevant beaches, the assumption of straight and parallel 
contours, and the use of a one dimensional model, appears justified. The 
"mean beach" surf information provided by the model will adequately 
represent the conditions over a standard 500 yard landing beach. If a wider 
beach needs to be characterized, the improved model would be run several 


times, entering the required information for each beach segment. 
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V. TEST AND EVALUATION OF THE IMPROVED MODEL 


A. WAVE HEIGHT AND BREAKER LOCATION 

The improved version of the SSSP was tested against one wave tank 
experiment and two sets of beach surf data to see if the model would choose 
reasonable, Gamma and B parameter values, based on the predictors discussed 
in Chapter 4. The predicted wave heights and surf characteristics are 
compared with the observed values. 

Case 5 (Battjes and Stive, 1985), discussed in the preceding chapter, 
was chosen as the wave tank experiment with which to test the finished 
model, This case was chosen because of its barred beach bathymetry. 
The model selected values of 0.41 and 0.963 for Gamma and B. The 
optimum values, from Table 4.1, are 0.45 and 1.30. The model wave height 
fit to the observed data points is still quite reasonable (Fig. 5.1). The 
solid line represents the model output and the circles are the observed data. 

Two Southern California beach data sets were used for further tests. 
Data collected at Santa Barbara's Leadbetter Beach during February, 1980 
and Torrey Pines during November, 1978 (Thornton and Guza, 1983), were 
used as the primary tests of the performance of SSSP. 

The first Santa Barbara test used data collected on 3 February, 1980. 
The average beach slope was 0.047 and the incident waves were 0.55 meters 
(rms) with a period of 14.3 seconds. The observed breaker height was 
measured at 0.70 meters, approximately 40 meters from the beach (dashed 


line). The wave angle was 7.8 degrees from normal. The SSSP predicted a 
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rms breaker height of 0.68 meters (solid line), 40.0 meters from the beach 
(Fig. 522) The model under predicted the observed heights less than 4 
percent, 

Similar resol are obtained for Leadbetter Beach data obtained on 4 and 
5 February, 1980. On 4 February, the incident wave height was 0.56 
meters with a period of 14.3 seconds. The wave angle was 9 degrees from 
normal, The observed breaker height was 0.67 meters, 50 meters from the 
beach. The model predicted 0.65 meters, 40 meters from the shore (Fig. 
ahs 

On 5 February, the incident wave was 0.45 meters with a period of 
12.8 seconds. The wave angle was 8.4 degrees from normal. Observed 
breakers were 0.60 meters, about 36 meters from the beach. SSSP predicted 
0.55 meter breakers, 33 meters from the shore (Fig. 5.4). 

Torrey Pines beach data acquired on 4 and 10 November, 1978 is 
compared with model results in Figures 5.5 and 5.6. The beach has a shallow 
slope of 0.022 which accounted for a wide surf zone of about 160 meters. 
On 4 November, the incident waves were about 0.35 meters with a period of 
14.3 seconds and normal incidence. The observed breaker height was 0.47 
meters about 110 meters from the shore. SSSP predicted 0.48 meter 
breakers at a distance of 106 meters from the beach. The model correctly 
forecast the wave height and surf zone width (Fig. 5.5). 

A similar result is seen for the November 10 Torrey Pines Case. The 
incident wave is 0.68 meters with a period of 15.9 seconds and normal 
incidence. Observed breakers were 0.72 meters, approximately 130 meters 
from the shoreline. The model predicted breaker heights of 0.75 meters, 140 


meters from the beach (Fig. 5.6). As for the November 4 case, the model 
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worked quite well for wave height prediction and the location of the breaker 
line. 

Based on the tests of the SSSP against various data sets, it appears that 
the new Person is an improvement over the original model. The old version 
consistently underpredicted the wave heights and was very inaccurate at 
estimating the location of the breaker line (Devendorf, 1985). The new 
version appears to be quite good at predicting the breaker height and gives 
better estimates of the breaker line. Table 5.1 is a comparison between the 
optimum values of the two breaker parameters and the values chosen by the 
model, The eaeaue for choosing Gamma and B is convenient, although 
improvements will need to be made in this area. 

B. LONGSHORE CURRENT DISTRIBUTION 

Longshore currents are the result of waves impinging at an oblique angle 
to the shoreline. All the wave tank Sesh. as well as the Torrey Pines 
cases, had normally incident waves. The Santa Barbara Leadbetter Beach 
data is the only comprehensive field data on longshore currents acquired to 
date. There exists no laboratory data for longshore currents generated by 
obliquely incident, random waves, The same three data sets (February 
3,4,5) were chosen as study cases because they had a variety of incident 
wave heights and wave angles. For the February 3 case (Fig. 5.7), the 
maximum observed current (dashed) was 0.48 m/sec at a distance of 22 


meters from the shoreline. The bed shear stress coefficient, C was set 


f’ 
equal to 0.009 for all cases. This value was derived experimentally by 
Thornton and Guza (1985). The SSSP prediction (solid) for longshore current 


was a maximum of 0.5 m/sec, 25 meters from the beach (solid line). 
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For the February 4 case (Fig. 5.8), the maximum observed current was 
.48 m/sec, 25 meters from shore; the model prediction was .44 m/sec at 
20 meters from the beach. The Feb 5 case (Fig. 5.9) shows an observed 
current maximum of .36 m/sec at 20 meters from the beach. The model 
predicted .36 m/sec at 20 meters from the shoreline, 

The model not only predicts the maximum current velocity reasonably 
well, but the current distribution within the surf zone shows good agreement 
with the data. The model seems to underpredict the current velocity 
seaward of the maximum current, but, in general, the model does quite a 
good job of longshore current prediction. 

The other Leadbetter Beach cases showed similar results. In all cases, 
the model gives reasonable maximum current locations and speeds as well as 
current distributions within the surf zones, These findings were also 
confirmed by Thornton and Guza (1985) in a-study in which the random wave 
model's wave height and longshore current predictions are rigorously tested 


against the Leadbetter Beach data. 
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VI. SUMMARY AND CONCLUSIONS 


A. SUMMARY 

In 1981, the U.S. Navy contracted to have a Sea, Swell and Surf 
Program (SSSP) developed to provide input to the planning efforts of Naval 
Oceanographers at sea and amphibious operational planners. The model was 
developed in FORTRAN and designed to run on a main frame computer (Wang 
and Chen, 1983). The completed SSSP was then translated into BASIC for 
use on the currently available shipboard micro-computers, the HP-9845B/275 
(Devendorf, 1985). The SSSP was modular in nature, with separate 
subroutines for offshore wind wave generation, surf zone calculations and 
nearshore circulation, 

The original SSSP had several problems which prevented its use by 
operational planners. The offshore wind wave generation module overbuilt 
the wave heights for high wind speeds. Some of the surf condition 
calculations were inaccurate and the nearshore circulation subroutine had 
stability problems. These problems were addressed by this thesis in an 
attempt to improve and upgrade the SSSP. 

The open water module shows satisfactory wind wave generation for wind 
speeds of 15 knots or less. For higher wind speeds, the predicted wave 
generation occurs too quickly, both spatially and temporally. Sensitivity 
analysis indicated that the Phillips linear term was responsible for the rapid 
development, but no improvements to the model physics were made. Several 
ad hoc "fixes" were recommended. This problem will require further 


research, 
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Two subroutines are called to calculate the wave transformation as the 
wave proceeds from offshore to onshore. Several numeric improvements were 
made which eliminated minor instabilities in the wave direction and height 
Proroutines. The original no flux boundary conditions were changed to cyclic 
boundary conditions in both’ subroutines, This change allows the wave 
heights and directions to more accurately reflect the conditions required by 
the model physics at the model boundaries. 

The most serious problem, the instability of the nearshore circulation 
module, is solved by substituting a one dimensional longshore current model 
for the original two dimensional scheme. Although there is a loss of two 
dimensional current resolution within the surf zone due to the assumptions of 
straight and parallel contours, stationary wave conditions, and shallow water 
initial waves, it was felt that a model that gives the important longshore 
current information quickly and accurately 1s an improvement over the 
original unstable, two-dimensional model. The new current model is a 
simple, longshore momentum flux balance equation. It retains the wave 
induced momentum flux and neglects the turbulent Reynold's stress. The new 
model is fast and accurate, requiring only a few calculations once the wave 
height field is generated. 

To take advantage of state of the art breaker height modeling, a random 
wave model was implemented for the calculations of the wave heights in the 
nearshore area. The random wave model assumes that wave heights, in the 
nearshore region, follow a Rayleigh distribution and that, even after the 
waves break, this distribution is retained.As a future project, the random 
wave model could be developed into a two dimensional model by applying the 


proper longshore numerical differencing. 
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The random wave model assumes that wave breaking can be described by 
a periodic bore. Two parameters are introduced by the model, which must 
be chosen based on the initial conditions. The Gamma term is part of the 
proportionality coefficient which modifies the initial Rayleigh distribution, and 
is found to be strongly correlated with the beach slope. The breaker 
parameter, B, describes the amount of the. wave face that is breaking and 
appears to be related to the surf parameter Eta. The criterion for choosing 
the two breaker parameters was developed empirically. 

The wave height and current models were next tested against wave tank 
and natural beach data sets. Although the model did not necessarily choose 
optimum values for the breaker parameters, the model is fairly robust, and 
slight deviations in the parameter values do not degrade the wave height 
calculations significantly. Good agreement was found between the observed 
wave heights and the model's predicted values. The longshore current, based 
on the model's predicted wave height field, were also in close agreement 
with observations. 

The calculation of Effective Surf, an amphibious planning guideline, was 
improved. The algorithm was implemented in accordance with the Joint Surf 
Manual (1976) with respect to its input requirements and its formatted 
output, 

A simple users manual, tutorial in nature, was developed for this version 
of the SSSP. It is short, and assumes some knowledge of the host computer 
and basic oceanographic terms. It is written within the guidelines set forth 
by the Navy for development of computer models for the Shipboard Numerical 


Aid Program (Brown, 1984). 
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B. CONCLUSIONS 

The improved Sea, Swell and Surf Program (SSSP) can provide much 
useful information to the Naval Qceanographer at sea or amphibious planner. 
It has several problems which have not been resolved, so its use_ is 
recommended with several caveats. The offshore wind wave generation 
module is felt to be accurate for low wind speeds only. Useful information 
on offshore wave growth can be obtained for input wind speeds of 15 knots or 
less. The results obtained for higher wind speed inputs must be examined 
carefully before use. The requirement for providing an initial wave energy 
spectrum is probably impractical for shipboard use. The option of providing a 
single initial wave height and period will be of more value to the operational 
user. 

The new random wave breaker height model appears to work well. The 
breaker heights and locations as well as the longshore current field are 
accurately predicted by the model. Due to its one-dimensional nature and 
the assumption of straight and parallel contours, longshore wave height and 
current resolution is restricted. However, the assumption that an amphibious 
landing beach is uniform in the longshore direction is a reasonable one over a 
wide range of operational areas. The small grid size of the model supports 
the use of a one dimensional model to compute an average set of charac- 
teristics for a "mean" beach. The speed and stability of the new model 
make it a viable tool for operational planners. Recommendations for further 
improvement of the SSSP include; 

1. Pursue the problem of the overbuilding of the offshore wind waves. 

2. Development of better predictors for the breaker parameters Gamma 


and B by including a wider data base. AS was seen in the 
optimized tests, correct choice of these parameters will provide a 
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very good fit to observations, even with very complicated 
bathymetry. 


3. Explore the feasibility of using the random wave model two 
dimensionally. This would offer the advantages of this more realistic 
approach to wave height modeling, while retaining the details that 
two dimensional spatial resolution provides. 

The improved SSSP offers the Naval Oceanographer at sea and amphibious 

operational planners a new and useful tool that, when used in conjunction 


with surf observations, weather information, and other available data, can be 


a great aid in the safe planning and execution of beach landing operations. 
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